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This paper presents a robust control design methodology based on the 
estimation of the first two order moments of the random variables and 
processes that describe the controlled response. Synthesis is performed by 
solving an multi-objective optimization problem where stability and perfor- 
mance requirements in time- and frequency- domains are integrated. The 
use of the first two order moments allows for the efficient estimation of the 
cost function thus for a faster synthesis algorithm. While reliability re- 
quirements are taken into account by using bounds to failure probabilities, 
requirements related to undesirable variability are implemented by quanti- 
fying the concentration of the random outcome about a deterministic target. 
The Hammersley Sequence Sampling and the First- and Second-Moment- 
Second-Order approximations are used to estimate the moments, whose 
accuracy and associated computational complexity are compared numeri- 
cally. Examples using output-feedback and full-state feedback with state 
estimation are used to demonstrate the ideas proposed. 
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I. Introduction 


Achieving balance between stability and performance in the presence of uncertainties is 
one of the fundamental challenges faced by control engineers. Trade-offs must be made to 
reach acceptable levels of stability and performance with adequate robustness to parameter 
uncertainty. These trade-offs are explicitly linked to the control engineer’s choice of uncer- 
tainty model as well as how that model is exploited in the synthesis process. Usually, the 
assumed uncertainty model has a profound impact on the performance robustness of the 
closed-loop system. 

Several uncertainty models, such as norm-bounded perturbations, interval analysis, fuzzy 
sets, and probabilistic methods 1-3 are typically used. The most commonly used robust 
control methods 4 are yU-synthesis and //-infinity. In these methods, uncertainty is modeled 
with norm-bounded complex perturbations of arbitrary structure about a nominal plant. 
This treatment is used primarily because it leads to a tractable set of sufficient conditions 
for robust stability, making the approach computationally efficient. These methods are based 
on the most pessimistic value of performance among all possible ones, usually referred to 
as “worst-case”. This worst-case performance is usually realized only by a single member 
of the uncertain model set and by a particular input signal. No information is provided 
regarding the likelihood that this worst-case will ever occur in practice. In addition, the 
intrinsic mathematical requirements of the approach usually lead to conservative models of 
uncertainty, over-conservative designs and complicated compensators. 

Probabilistic uncertainty not only defines a set of plants where the actual dynamic sys- 
tem is assumed to reside but also associates a weight, the value of the probability density 
function, to each member of the set. In contrast to conventional robust control methods, this 
“additional dimension” allows the pursuit of robustly optimal solutions in the probabilistic 
sense. For instance, reliability-based design searches for solutions that minimize the proba- 
bility of violating design requirements prescribed in terms of inequality constraints. Hence, 
reliability-based control design searches for the compensator that places as much probability 
as possible within the region where the design requirements are satisfied. Notice that this 
allows for the search of the compensator with the best robustness for a given control struc- 
ture, e.g. the most robust PID controller. This can be achieved even though the violation of 
some the design requirements for some of the plants in the uncertainty set is unavoidable. 

Synthesis approaches based on random searches 5, 6 ’ 14 and stochastic gradient algorithms 8-10 
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have been applied to probabilistic robust control. In these studies, random sampling is the 
primary tool for assessing and pursuing acceptable levels of robustness in the control solution. 
On the other hand, asymptotic approximations 11,12 for the estimation of failure probabili- 
ties have only been used as a control analysis tool. Reference, 13 foundation of this study, 
integrates and extends some of these tools. 

This paper is organized as follows. Section II presents basic concepts related to control 
and probabilistic uncertainty. Section III introduces reliability metrics for random variables 
and processes. Mean and variance based bounds to the reliability metrics are also derived 
therein. Robustness-based metrics for random variables and processes are introduced in 
Section IV. Section V presents the numerical methods used to estimate the above mentioned 
metrics. The control synthesis procedure is presented in Section VI, where specifics of 
both the reliability and the robustness-based formulations are examined. Two examples are 
presented in Section VII, where a satellite’s attitude control problem and the disturbance 
rejection in a flexible beam are used to demonstrate the method. Finally, some conclusions 
are stated in Section VIII. 


II. System Dynamics 

Let p be a vector of random variables used to model the uncertain parameters of the 
system. In this study, p is prescribed a priori by the joint probability density function 
(PDF) / P (p) or equivalently by the cumulative distribution function (CDF) F p (p) a . The set 
of values that p could take, called the support of p, will be denoted as A p . 

Consider the probabilistic model AI(p) of a Linear Time Invariant (LTI) system, where 
the dependence of the model on the uncertain parameters could be non-linear. The reader 
must notice however, that the developments presented herein do not require the system to 
be LTI. The propagation of A p through A4 leads to a set of uncertain plant models in which 
the physical system is assumed to reside. The probability distribution of a plant within this 
set is fully determined by -M(p) and / p ( p). In a transfer function representation, we will 
refer to the uncertain plant as G(p) and to the compensator as iC(k), where k is the vector 
of design parameters to be determined. Alternatively, a state space realization of Ad(p) 

a In these expressions, the subscript refers to the symbol used for the random variable while the value in 
parenthesis refers to a particular realization. 
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leads to 


( 1 ) 

( 2 ) 


x = A(p)x + B(p)u + F(p)z 
y = C(p)x + D(p)u + E(p)v 

where x is the state, u is the control, z is process noise, y is the system output and v is 
sensor noise. The noise signals are commonly modeled as delta correlated Gaussian white 
noises satisfying E[z] = 0 and E[z(f)z 2 (t + «)] = S<5(/c), where z = [z T ,v T ] T and S is a 
constant spectral density matrix. In what follows, the explicit dependence of the matrices 
in Equations (1-2) on p is omitted while D is assumed to be zero. 

Important properties, typically used in control design, such as pole placement and the 
Separation Principle, do not hold due to the offset between the deterministic mathematical 
model and the actual dynamic system. The effects of parametric uncertainty on the Separa- 
tion Principle are considered next. For the full-state feedback law u = — Gx and a full-order 
observer with gain L based on the expected plant E[At(p)] (any other deterministic plant 
such as _Ad(E[p]) could be used instead), the closed-loop dynamics for full-state feedback is 
given by 


x = Ax + Bz 

(3) 

y = Cx + Ez 

(4) 



A - BG 

BG 

A = 

A — E[A] + (E[B] — B)G 
+L(E[C] - C) 

(B — E[B])G + E[A] — LE[C] 


B 


F 

0 

F 

-LE 


where x is the estimation of x, x = [x T ,e T ] T is the augmented state vector, e = x — x 
is the estimation error, C = [C 2 0 2 ] T and E = [0 2 |E 2 ] 2 . The vector k is formed by the 
feedback gain G and the observer gain L. Notice that the Separation Principle holds, i.e., 
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A is upper triangular, if the deterministic plant used to generate the observer matches the 
physical dynamic system. Uncertainty in the plant makes the Separation Principle invalid. 
In addition, the random closed-loop poles do not occur at the locations selected for the 
full-state feedback, i.e. poles of the subsystem, nor at the locations for the full-order 
observer, i.e. poles of the A 2 ,2 subsystem. 

III. Reliability-Based Metrics 

The propagation of a fixed set of parameters of the plant through conventional control 
analysis tools leads to set of scalar quantities, e.g. closed loop poles, and a set of functions, 
e.g. step responses, and Bode plots. The propagation of probabilistic uncertainty through the 
same tools leads to random variables, e.g. random closed-loop poles, and random processes, 
e.g. the step responses become random processes parameterized by time, and the Bode plots 
become random processes parameterized by frequency. In this section we first introduce 
reliability metrics for random variables and processes. These metrics will be used to quantify 
the violation of the design requirements. Specific realizations corresponding to stability, time, 
and frequency requirements are then provided. In general, we will use x and x(h) to denote 
a random variable and a random process dependent on p through the plant model. For the 
random process x(h), h refers to an arbitrary variable such as time or frequency. 

A. Random Variables 

We first introduce the concept of probability of failure. Let x(p) be the random variable of 
interest. Let x > x be a design requirement. The event x < x will be referred to as failure. 
The corresponding failure set is given by T = {x \ x G (— 00 , 2 ;]}, where the failure boundary 
x is a deterministic quantity prescribed in advance. The admissible domain, namely A = 
{ 2 ; j x G (x, 00 )}, is the complement of the failure domain. The same type of discrimination 
can be done in the parameter space p by using x(p). The function g( p, x) = x(p) —x, called 
the limit state function , divides the parameter space in two parts, the domain leading to 
A, i.e. g(p,x) > 0, and the domain leading to A, i.e. g(p,x) < 0. Hence, T results from 
mapping the set {p G A p | g(p,x) < 0} through x(p). In this case, the probability of failure 
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Pf is given by 


P f = P[x < x] = f f x [x)dx = f f P (p)dp (5) 

J x<x J g< 0 

Similar expressions can be derived if the design requirement is x < x. These expressions 
describe reliability metrics for the random variable x when a single constraint is present, i.e., 
x > x or x < x. A reliability metric for the random variable x having both constraints can 
be easily formed 

r x (x, x) = r x (x) + r x (x) (6) 

where 

L x (x) = P[x < x] = F x (x) (7) 

r x {x) = P[x > x] = 1 - F x (x) (8) 

Notice that rjx) is equivalent to Pf in Equation (5). We will refer to x and x as the 
boundaries of the failure domain F = {x \ x G (— oo,;r] U (x, cxd)}. Notice that the under-bar 
and the over-bar refer to the bound from below and the bound from above of the admissible 
domain A = {x j x G [x,x)}. This convention will be used for the remainder of the paper. 
Notice that the mapping of the corresponding limit state function through x(p) leads to the 
failure boundary (s). Hence, there is a direct correspondence between T and g. 

B. Random Processes 

The random process x(h) can be considered as the parameterization of the random variable 
x by the deterministic variable h. In this paper h G [0, oo] is assumed. The random process 
x(h) is specified by the set of CDFs 15 F x ^(x,h). For instance, the system output y(t) is 
prescribed by F y ^)(y,t)- The evaluation of the process at a particular h value, say hi, leads 
to the random variable x(hi) whose CDF is given by F x (x) = F x ^(x,hi). In general, the 
support and the percentiles of x(h) vary with h. 

Let x(h) > x(h) for h G [hi, I 12 ] and x(h) < x(h) for h G [h 3 , h^] be design requirements 
for the random process x(h). In this paper, reliability metrics for processes are formulated 
by extending the ideas presented above. This is attained by integrating the reliability metric 
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in Equation (6) for the random variable x{hi) in the h — interval of interest. In this context, 
a reliability metric for x(h) is cast as 

r x ( h ) (x{h),x(h)) = r x{h) (x(h)) +r x{h) (x(h)) (9) 

where 

r x{h) (x(h)) = - — ] —r~[ P [x(h) < x{h)]dh = - — 1 —l F x{h) (x(h), h)dh (10) 

ft 2 — ft i Jhi ~ m Jhi 

r x(h)ix{h)) = - — ^-r-f P [x(h) > x(h)]dh = - — ^-r- [ 1 - F x{h) (x(h),h)dh (11) 

h 4 — Jh 3 /? '4 — ft 3 .7/1,3 

are the costs of violating the lower and upper constraints respectively. These constraints, 
namely x(h) and x(h), will also be referred to as failure boundary functions. Notice that the 
failure domain 


T = { 1J {(x,h)\x < x(h)} 
he[hi,h 2 ] 


U 


U {{x,h)\x>x{h)} 
he[h 3 ,h 4 ] 


is delimited by the failure boundaries. The reader shall realize that Equation (9) is a natural 
extension of Equation (6). If the process is contained within the set A the reliability metric 
r x (h) is zero, meaning that the inequality constraints are satisfied for all parameter values in 

A P - 


C. Realizations 

1. Robust Stability 

A LTI system is robustly stable if all its poles are in the open left half of the complex plane 
for all possible values of the random parameters. A reliability assessment of stability is given 
by 


P 


U(»M > o) 

. 2=1 


6 


where Sj with i = 1, 2, ... u is a random pole, SR[-] is the real part operator and e is the 
resulting probability of instability. Robust stability is attained if e = 0. Stability can also 
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be cast via 


A = max{9ft[si], 9ft [S 2 ], . . . , 9ft[s.y]} (12) 

In terms of A, the probability of instability is given by r_\(0). Robust stability is attained 
if Fa( 0) = 0. Several comments are now pertinent. Reaching robust stability may not be 
feasible for the given support A p (even if it is bounded) and the assumed control structure 
K( k). Notice also that the acceptance of a small non-zero probability of instability could be 
desirable from the performance point of view. For instance, allowing the right low-probability 
tail of /a (A) to lie on the open right half of the complex plane may yield a significant 
enhancement in the performance of the plants associated with the high probability portions 
of the PDF. Rather than advocating for the acceptance of the risk that this practice implies, 
we would like to highlight that the trade-off between robustness and performance can be 
studied by allowing small values of e. 

2. Time- Domain 

Quite frequently performance requirements are prescribed in terms of time-domain specifi- 
cations. The propagation of / P (p) through the system dynamics leads to random processes 
for the time responses. Denote by x(t) an arbitrary random process with CDF F x u^(x,t). 
Such process is parameterized by p, time t, and the compensator design variable k. The 
dependence of x(t) on k has been omitted for the sake of simplifying the notation. Reliability 
metrics for relevant processes can be cast using Equation (9). For instance, while settling 
time and overshoot requirements are integrated using r y ^(y(t),y(t)), the control saturation 
requirement \u\ < u max leads to r u ^(-u max ,u max ). 

A reliability metric for assessing the effects of noise on the uncertain plant is formulated 
next. The state covariance matrix, defined as Q (t) = E[x(f)x T (f)], is given by the solution 
to the covariance equation 


Q = AQ + QA t + BSB t 


(13) 


subject to Q(0) = Qo- The output covariance, defined as E[y(f)y T (f)], reaches the steady- 
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state Root Mean Square (RMS) value 


y rms = lim diag 


CQ(t)C 7 


1/2 


(14) 


Notice that uncertainty in p makes y rms a random vector. If y rms is a component of y rms , 
a reliability metric that penalizes the violation y rms > y rms is given by r yrms (y rms ). 


3. Frequency-Domain 

The propagation of / p ( p) through the system dynamics onto the frequency domain leads to 
random processes of the form x(co), fully specified by F x ^)(x, u). Here, x(u) is any real fre- 
quency dependent metric of the feedback loop, e.g. Bode magnitude. This random process 
is parameterized by p, frequency u, and the compensator design variable k. A reliability 
metric for x(u) is r x ( u ,)(x( u>),x(u>)). For instance, conventional control requirements 16 for dis- 
turbance rejection, noise attenuation and reference tracking can be cast in terms of the loop 
transfer function q( u) = \GK\. Low frequency requirements can be cast using r qH (q{u)) 
with q(u>) = 1 and high frequency requirements with r q ( u )(q( u)) for which q( u) has a proper 
roll off. 


D. Reliability Bounds 

The following lemma defines bounds for the reliability metrics in terms of the first two order 
moments. 

Lemma 1 . Let E[-] and V[-] denote the expected value and variance operators. 


L x {x) < 

r x {x ) < 


V[x] 

(E[x] — x) 2 
V[x] 

(x — E[a:]) 2 


if x < E[x] 
ifx> E[x] 


(15) 

(16) 


Proof. If e > 0 


V[x] > f (x — E[a;]) 2 f x (x)dx > e 2 P [|a; — E[x]| > e] > e 2 P [x < E[x] — e] 

J \x—E[x]\>e 

We obtain Equation (15) using x = E[x] — e. Equation (16) is derived following the same 
lines. □ 
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Notice that these bounds apply to any PDF of x. Since the bounds given in Equations 
(15-16) are solely dependent on the first two order moments they can be estimated more 
efficiently than exact failure probabilities. Notice however that the use of bounds instead of 
reliability metric introduces conservatism into the solution. 


where 


b x (x. 

<x) =b x (x) 

+ b x (x) 

(17) 


f v[*i 

if x < E[x] 
otherwise 


b x {x) = < 

1 (E[a;]-s) 2 

OO 

(18) 

| 

f V[*] 

if x > E[x] 
otherwise 


b x (x) = < 

1 (®-E[q) 2 

[ 00 

(19) 


As before, the under-bar and over-bar in b refer to the way in which the failure event is 
defined. Notice that r x < b x , r x (x) < b x (x), and f x (x) < b x (x). The same idea, extended to 
random processes, leads to 


(x(h),x(h)) = b, m (x(-)) + b xm (x(h)) 


( 20 ) 


where 


a I b x (x(h))dh if x(h) < E[x(/i)] Mh e [hi, h 2 


b x{h) (x(h)) = 


oo 


otherwise 


( 21 ) 


b x (h){x(h)) = 


a I b x (x(h))dh if x(h) > E [x(h)} Mh e [h 3 , h 4 \ 


oo 


otherwise 


( 22 ) 


As before, r x{h) < b xW , r x{h) (x(h)) < b x{h) (x(h)) and r x{h) (x(h)) < b x{h) (x(h)). 


IV. Robustness-Based Metrics 

Some performance requirements might not be properly captured in a conventional relia- 
bility formulation since they focus on the bulk portion of the PDF. Throughout this paper 
the term robustness is used to describe the design characteristic that evaluates the perfor- 
mance degradation from an ideal deterministic behavior caused by uncertainty. Robustness 
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metrics, that quantify such a characteristic, are presented next. For random variables, the 
index 

r x (x) = f (f -x) 2 f x (£,)d£ = V[a] + E[x](E[x] - 2x) + x 2 (23) 

J A x 

is a measure of the concentration of f x (x) about the deterministic target value x. In the 
ideal case when f x (x) = S(x — x), we obtain r x (x) = 0. For random processes, the index 

1 r K 

r x (h)(x(-)) = t~ / r x (x(h))dh (24) 

h 5 Jh 5 

is a measure of the concentration of the process x(h) about the deterministic target function 
x(h) in h E The reader must notice that the evaluation of the above expressions 

only requires of the first two order moments of the process. 

Robustness-based metrics parallel to the ones provided in Section C can easily be posed. 
For instance, if y rms is the RMS steady state value of an error signal, the index T yrms ( 0) 
quantifies the offset between the target behavior y rms = 0 and the random variable y rms . 
Likewise, the metric r u ( t )(0) quantifies the offset between the random process u(t) and the 
target u — 0, for which no actuation is required. 

V. Numerical Estimation 

Since the estimation of the metrics for random processes are a natural extension of the 
methods for random variables, only the estimation of the statistics of x is addressed here. 
The extension to processes is as follows. For the random process x(h), create a uniform 
sample of e points hi, h 2, . . ., h e in the /i-domain. We will refer to these samples as the e 
h-samples. Statistics for the resulting e random variables x t = x{hi), i = 1, 2 ... e are then 
used to estimate the pertinent integrals, i.e. Equations (10-11), (21-22) and (24). 

A. Mean and Variances 

The reliability bounds and the robustness metrics introduced above depend exclusively on 
mean and variances. In this paper, Sampling and the First and Second Moment Second 
Order approximations are used for their estimation. These methods are briefly introduced 
next. 
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1. Sampling 

Unbiased estimators for the mean and variance of the random variable x are 

1 n 

E[x] « — ^ Xi (25) 

n i = i 

n 

v[x] « - E[x]) 2 (26) 

n — 1 

1=1 

where Xi = x(pj) is the ith sample. The p, sample of / P (p) can be generated by any sampling 
technique. 

HSS generates representative deterministic samples of / P (p). The error of approximating 
an integral by a finite number of samples of the integrand, depends on the uniformity of 
the points used to generate the samples rather than on their randomness. This fact has 
motivated the development of deterministic sampling techniques such as HSS. While con- 
ventional Monte Carlo Sampling (MCS) is based on the generation of random points on the 
unit hypercube, HSS is based on the generation of an evenly distributed set points. The 
reader can refer to 13 for details on HSS implementation. 

HSS requires far fewer samples 18 than MCS 5 ~ 7,14 for a given confidence level. Improve- 
ments in the convergence rate of the estimated first two order moments by a factor of three 
to one hundred 19 have been reported. In addition, the estimated value of the failure proba- 
bility based on HSS is deterministic. In contrast, MCS leads to random estimates unless an 
infinite number of samples is used. This is especially noticeable if n is small. The random 
character of the estimation can only be mitigated by increasing the number of samples, which 
incidentally increases the computational demands of algorithms based on MCS. Therefore, 
HSS not only leads to more accurate estimations than MCS for a given number of samples 
but also eliminates the random character of the results. 

2. First and Second Moment Second Order Approximations (FSMSO) 

These approximations result from calculating the first and second moment of a second order 
Taylor expansion of x(p) about E[p]. Let p G M m , and p, be the ith component of p. If the 
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components of p are independent random variables, the resulting approximations are 


i AL f ) 2 r 

E M ~ iTTT V [P«] 


v[ x ] « J2 


i= 1 
m m 


br a p 

dx 

^P(i) 


(0 


V[ P( i)] 


dx 

<9pp) 


d 2 x 


Dp 2 


T [P(d] + 7 


d 2 x 


W, 


4 \ dp 2 


« 


SSKsSs) V| "' |VM 


Y 


*=i 


( 27 ) 


(F[p (i) ] - V[p (t| ] 2 ) 
(28) 


where the functions and derivatives are evaluated at E[p], T[-] is the third central moment 
operator and F [-] is the fourth. 

In general, reliability metrics cannot be evaluated exactly since they involve the evaluation 
of complicated integrals, usually multi-dimensional, over complex domains. The estimation 
of failure probabilities, basic component of the reliability metrics, can be done using sampling 
or asymptotic approximations. The reader can refer to 13 for details on their estimation. 


VI. Control Synthesis 


A. Reliability-based 

The formulation of the control design problem from a reliability perspective is as follows. 
For a given plant model, compensator structure, uncertainty model, and a set of design 
requirements prescribed via inequality constraints; one would like to find the compensator 
parameters for which the resulting probability of violating the design requirements is mini- 
mized. Notice that this refers to the excursion of the outcomes into the failure domains. 

1. Robustness Considerations 

The reliability metrics in Equations (6-9) are usually applied using a fixed failure set T . 
In this form, a reliability analysis cannot assess the system’s performance in the regions 
where the design requirements are satisfied, i.e. the intersection of the admissible domains 
associated with all the design requirements. Since the portion of the random outcome lying 
in the admissible domain A might end up being substantially larger than the portion lying 
in the failure domain F, a reliability-based approach with fixed failure boundaries does not 
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have control over the bulk portion of the PDF, which is the portion that dictates the most 
likely performance. 

We now introduce the concept of a shapeable failure set T with an example. Let x(k) 
be the stationary RMS value of an error signal. One would like to find k such that x is as 
close as possible to zero. Uncertainty in the plant makes x a random variable. Let x be 
the failure boundary associated with the design requirement x < x, i.e. a fixed failure set is 
J- = {x | x G [x, oo)}. The minimization of r x (x) leads to reliability optimal compensators. 
Suppose there exist multiple designs leading to r x (x) = 0. These designs however differ in 
how well the resulting PDF of x spreads over the admissible domain A = {x \ x G [0, x)}. 
The concentration of f x (x) about zero is an indicator of the robust performance. Say, for 
example, that kx leads to r x (x/ 2) = 0 and k 2 leads to r x (x) = 0. Since neither of these 
two designs violate the design requirement, a reliability analysis cannot establish that the 
compensator with parameters ki has a better robust performance than the one which uses 
k 2 . 

By minimizing the reliability metrics and simultaneously shrinking the admissible do- 
main, the whole random variable/process can be concentrated about regions with an im- 
proved system performance. This is attained by parameterizing the failure boundaries of 
T as well as a penalizing function 7 X with an additional design variable, namely e. The 
basic idea is to solve an optimization problem for the design variable d = [k, e] such that a 
twofold objective is pursued: a reliability metric is minimized while the size of A is reduced. 
For the RMS example above, the minimization of J = r x (e) + j x where 7 X = e, d = [k, e] 
and e e [0, x], leads to the desired solution. This setting implies T = {a: | x G [e, oo)} and 
A = {x | x G [0, e)}. See Figure 1. Notice that the value of J for ki is less than the one for 
k 2 if e G [x/2, x). In general, we will refer to the augmented reliability metric as the sum of a 
reliability metric from Section III and a penalizing term. Augmented reliability metrics for 
the random variable x and the random process x(h) take the form r x (x(e),x(e)) + 7 x (e) and 
r x{h)(%(h,e),x(h,e)) +'y x (h)( e ) The penalizing functions 7 x (e) and 7 x (h)( e ) must be propor- 
tional to the size of the admissible domain A. In addition, they must be built such that the 
minimization of the augmented metric does not lead to unacceptable solutions, e.g. r x = 1 
and 7 X = 0. If r x < e is required, use a monotonically increasing function satisfying 7 e [0, e]. 
In the RMS example above, this is attained by minimizing an augmented reliability metric 
with 7 x (e) = ee/x for e e [0, x\. 
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Figure 1. Robust performance concepts and shapeable failure domains 


B. Robustness-based 

The formulation of the control design problem from a robustness perspective is as follows. 
For a given plant model, compensator structure, and uncertainty model; one would like to 
find the compensator parameters for which the resulting random outcome is as concentrated 
as possible to a target deterministic behavior. While the deterministic target for a random 
variable is a scalar, the target for a random process is a function, e.g., u{t) = 0 is the target 
function for the control. 

C. Synthesis Procedure 

A step-by-step procedure to control synthesis is presented next. 

1. Determine the plant model and the control structure. First principles and classical 
deterministic approaches to compensator design can be used. Identify the set of pa- 
rameters that have a strong impact on the plant model. Use sensitivity information 
and engineering judgment to select the set of uncertain parameters p. At this stage, 
the parametric plant model, i.e. G(p), and the control structure, i.e. K{ k), must be 
fully determined. 
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2. Generate the probabilistic parameter model / P (p). Use engineering judgment and 
experimental data if available. 

3. Cast the design requirements, either reliability or robustness based, in terms of the 
metrics introduced above. Use Equations (6,9) for the reliability metrics and Equa- 
tions (23,24) for the robustness metrics. Recall that while each reliability requirement 
requires setting a failure domain JE, each robustness requirement requires setting a 
target behavior. Use these metrics to compose the cost vector c, which is the vector of 
objectives for multi-objective optimization. 

4. Let d be the design variable. For robustness-based metrics and reliability-based metrics 
with fixed failure domains, d = k. Reliability-metrics with shapeable failure domains 
lead to d = [k T ,e T ] T . Build a penalizing function 7(e) for these terms. Llpdate the 
components of the cost vector c by adding the penalizing functions and parameterizing 
the failure boundaries. 

5. Solve the single objective optimization problem 

d* = argmin |c T w} (29) 

where w is composed of non-negative weights. Each cost function evaluation used in 
the search for the optimal design d* requires a probabilistic analysis. This analysis is 
done by calculating the metrics contained in c. Either sampling or FSMSO can be used 
to estimate the mean and variance based metrics introduced above. This task requires 
forming the closed-loop Equations (3-4) and performing typical control studies such as 
hireling closed loop poles, time responses and Bode plots. 

The cost vector c can be formed by combinations of reliability-based metrics, reliability 
bounds and/or robustness-based metrics. The corresponding implications are explored in 
the examples. 

D. Optimization under Uncertainty 

Let J(p,d) = c 7 w denote the cost function of the optimization problem. If reliability 
metrics are used in in c, the cost function might not only have plateaus, i.e., there could 
exist a design d and a non- zero perturbation 5 such that J(p,d) = J(p,d + <5), but might 
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also have a discontinuous gradient. The use of sampling in the estimation of probabilities 
makes the cost function piecewise constant. Let J e ( p, d) be an estimation of the actual cost 
J(p,d). For any design d and regardless of the number of samples, there always exists a 
perturbation § such that J e (p,d) = J e (p,d + 5). This situation is aggravated, i.e. bigger 
perturbations can be found, when a smaller number of samples is used or when Pf is close 
to zero or one. 

Mean and variances can be estimated via HSS or FSMSO. Their estimation via sampling 
does not exhibit the numerical problems mentioned above. Among the sampling techniques, 
it was found that HSS suites very well the need for an accurate and efficient estimation. 18 
Estimation via FSMSO is considerably faster than sampling but less accurate in general. 
This is because, the Taylor expansion might become a poor approximation of the actual 
function when / p ( p) is not concentrated about E[p], The methods used for the estimation 
of J e must be taken into account when selecting an optimization algorithm. In this paper, we 
use the Nelder Mead Simplex algorithm, which is a local non-gradient based search method. 

VII. Numerical Examples 

The synthesis procedure of Section C is applied herein to two examples. A textbook 
satellite attitude control problem is considered first. Then, the active control of a flexible 
beam subject to disturbances is considered. The following parameter values are assumed in 
both problems. If p G M m , and sampling is used, a coarse probabilistic analysis is performed 
using n = 75 m p-samples and e = 90 h— samples. If the resulting assessment is satisfactory, a 
new analysis is performed using n = 500m p-samples and e = 180 h— samples. This two-step 
analysis is carried out to adjust the fidelity of the estimation according to the characteristics 
of the compensator. For the sake of comparison, the examples also present the analysis of 
deterministic versions of the problems for which E[p] is used as parameter values. Such 
problems and their solutions will be referred to as the nominal ones. 

A. Satellite Attitude Control 

Accurate satellite pointing in the presence of large thermal gradients and mass losses for 
uncertain initial conditions is desired. A simple rotational model of two bodies connected 
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with a flexible boom leads to 


J\0\ + b(9\ — 62) + k{9\ — 62) — u 
J2O2 + b(d 2 — Oi) + k{9 2 — 0±) = 0 


where 0\ and 0 2 are the deflection angles, J\ and J2 are moments of inertia, k is the equivalent 
stiffness, b = a^/k / 10 is the equivalent damping coefficient and u is an applied torque. The 
variable a is used to model the changes in damping caused by thermal variations. Assume 
that J2 = 0.1 since mass losses only affect J\. The non-collocated sensor-actuator pair 
resulting from using y = 6 2 leads to the SISO system 


G(p) 


k + bs 

did 2 s 4 + b(J\ + t/ 2 )s 3 + (Ji + J2)ks 2 


(30) 


Variations in the operating conditions and ignorance on the initial conditions are modeled 
using p = [Ji, a, k, 6 W , 6 W , O 20 , d 20 ] T , where 6 W = 0i(O) and O 20 = d 2 ( 0). The following 
output-feedback control structure is assumed 


K{ k) 


k\ + /c 2 s + k^s 2 + /C4S 3 
k 5 + k 6 s + k 7 s 2 + k$s 3 


(31) 


The joint PDF that describes the uncertainty in p is given by the independent random 
variables listed in Table 1, where U and B refer to Uniform and Beta distributions. Notice 
that the Beta distribution has four independent parameters, two of them are the conventional 
arguments and the other two are the support bounds. 

1. Nominal Compensator 

A baseline compensator for the nominal plant is designed by standard pole placement tech- 
niques such that the resulting closed-loop system is robustly stable for the uncertainty set 
defined above. This was achieved by pursuing large stability margins. This practice results 
in a compensator with parameters kx = 10 6 [0. 0108, —0.3271, 0.1192, 0.0092,1.8835,2.1305, 
2.2276, 0.9308] t . The analysis of the nominal compensator using / P (p) indicates that the 
closed-loop system is robustly stable as intended, i.e. rA(0) = 0, but the time responses are 
unsatisfactory. The CDF of A as well as the time evolutions of the output and the control 
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Table 1. Uncertainty model. 


•/, 

A* = [0.8, 1] 

= 0.8,1) 

a 

A a = [0.03,0.2] 

f a (a) = B( 0.3, 0.2) 

k 

A k = [0.09, 0.4] 

fk{k) = B{ 5,5) 

#10 

A» 1( , = |-*/2,w/2] 

f 01o (6 w ) = B( 5.2, 5.2) 

#10 

A *„ = [-15,15] 

4 O (0 1O ) = B( 2.5, 2.5) 

#20 

Ae 20 = [— tt/2, tt/2] 

/ 02O (#2o) = 5(5.2, 5.2) 

#20 

4 2c1 = [-15.15] 

4,(020) = -0(2.5, 2.5) 


signals are shown in Figures 2-4. The sudden variation in the slope of the CDF of A in Fig- 
ure 2 is the result of a change in the closed-loop pole that determines A. The considerable 
disparity between A(E[p]) and E[A(p)] shows that the nominal problem is not a meaningful 
representative of the probabilistic behavior. Figures 3 and 4 show the time evolution of the 
random signals by indicating the 1, 10, 20, 30, 40, 50, 60, 70, 80, 90 and 99 percentiles. In 
Figures 3 and 4, the percentiles and the nominal functions are shown. Dotted lines are used 
to indicate the failure boundaries, specified in the next section. It is interesting to see how 
the PDFs expand, e.g. Figure 4 at 2.5 and 8 seconds, and contract, e.g. Figure 4 at 4 and 
16 seconds, in a oscillatory manner. This information can be used to determine the time 
periods when the effects of uncertainty are more noticeable. 

2. Mean and Variance based compensator 

If only means and variances are used, reliability and robustness metrics can be combined to 
form c = [6a (0), T y{t) ( 1), r u(t )(0), b q{u) ( 1), b q ^)(0. 75 / lj)] t . Ranges provided before will be 
used with the exception of uj G [5, 10 2 ] for b. The cost vector can be estimated via sampling 
or the FSMSM. 

HSS leads to k 3 = 10 5 [1.1207, 0.2307, 1.5235, 0.1713, 0.001, 0.0692, 1.2294, 1.069] T and c = 
[8.6 x 10 -5 , 0.0057, 0.0129, 1.94 x 10 -4 , 6.189 x 10~ 5 ] T for which the weighting vector w = 
[500, 1,4, 1, 1] T was used. A probabilistic analysis of this compensator leads to Figures 5, 6 
and 7. Figure 5 shows that the optimal compensator makes A insensitive to uncertainty, i.e. 
the CDF resembles a step function which would result if the system is deterministic. The 
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*(E[p]) E[X.(p)] 



Figure 2. CDF of A for the nominal compensator. 


reader must notice that the bound on the stability condition tends to reduce the variability 
of A. This artificially reduces the design space, i.e. there might exist solutions with a better 
performance in which large variances in A occur. This is a consequence of the form of the 
bound. This problem has been studied in reference 13 using reliability metrics exclusively. 
We will refer to this compensator as the reliability compensator. The comparison between 
the time responses for the reliability compensator and Figures 6 and 7 show substantial 
differences between the two solutions. While the control for the reliability based compensator 
intends to keep the process within the strip |u| < 0.5, the mean and variance based solution 
intends to concentrate the random process about the target function u — 0. It can be seen 
that this is achieved by inducing oscillations about the target function. The reader should 
also notice that excursions beyond the failure boundaries in a reliability formulation, are 
not penalized according to the severity of the violation. This is in sharp contrast with the 
mean and variance formulation. The CPU time for the finer assessment, i.e. the one using 
n = 500m p samples, is 77s. For a given number of samples, the estimation of means and 
variances is in general more accurate than the estimation of failure probabilities. Small 
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Figure 3. y(t) for the nominal compensator. A zoom is shown below. 


sample sets however, result in larger estimation errors of the moments, which could lead to 
the selection of the wrong conditions in Equations (18-19) and (21-22). 

Next, the FSMSO method is used. The first and second order derivatives for all the 
metrics of interest, i.e. closed-loop poles, output, control and loop transfer function, were 
derived analytically. Some of the required expressions are as follows 


X 


A 

0 

0 


X 


B 

<9x 

= 

<9A 

A 

0 


<9x 

+ 

<9B 

<9 2 x 


d 2 A 

2<9A 

A 


<9 2 i 


<9 2 B 










y 


COO 


X 

dy 

= 

dC C 0 


<9x 

d 2 y 


d 2 C 2dC C 


Qd 

to 

X' 


(32) 


(33) 


where A, B, C, and D were given in Equations (3) and (4), x = [xJ,x 7 ] T , <9[-] indicates a 
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Figure 4. u(t) for the nominal compensator. 


derivative with respect to p(j) excluding initial conditions, and 


A = 


A c 

— B c C 

BC c 

A - BD c C 


<9 A = 


0 -B c dC 

dBC c dA — BD c <9C — <9BD c C 


d 2 A = 


0 


—B r d 2 C 


d 2 BC r d 2 A - BB r d 2 C - 2dBB r dC - d 2 BD r C 


B = 


B c BD C 


, 5B = 


0 <9BD r 


, d 2 B = 


0 d 2 BD f 
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Figure 5. CDF of A for the mean and variance based compensator. 


c = 


0 c 


, dC = 


0 dC 


, d 2 c 


0 d 2 C 


In these expressions, the subscript c refers to the state space representation of the com- 
pensator K while the matrices with no subscript refer to the open-loop plant. The time 
evolution of the sensitivities is calculated by solving the state space model in Equations 
(32-33). Sensitivities with respect to initial conditions can be analytically calculated via the 
matrix exponential. Such developments as well as the ones for the other sensitivities are 
omitted due to space limitations. In general, the analyses based on FSMSO were inaccurate 
when compared with sampling. Even though the CPU time per analysis was reduced to 
2 seconds, large errors in the estimation precluded its use for synthesis. Table 2 shows a 
comparison between the assessments resulting from HSS and FSMSO as the variance of the 
uncertain parameters is increased. The same PDFs for the input uncertainty of Table 1 are 
used while the means are kept constant. The right most column shows the average relative 
error in the components of c. For the Beta distributions, the increase of the variance was 
attained by enlarging the support. The terms resulting from the crossed derivatives were 
neglected. It can be observed that the accuracy of the FSMSO estimation rapidly decreases 
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Figure 6. y(t) for the mean and variance based compensator. A zoom is shown below. 


as / P (p) is less concentrated about its mean. This trait is obvious since large excursions 
from E[p] might considerably degrade the accuracy of the Taylor approximation. 

The use of reliability bounds introduces unnecessary conservatism into the problem. To 
show the reduction in the design space caused by using & A (0), we have solved this prob- 
lem using Fa(0) instead. This practice results in the compensator k 4 = 10 5 [0.9082, 0.3590, 
1.1773, 0.2372, 0.001, 0.0616, 1.0432, 1.196] T for which c = [1.86 x 10” 5 , 0.0052, 0.0121, 1.845 x 
10 -4 , 6.10 x 10 5 ] T . Figure 8 shows that A a is about 10 times larger than the one in Figure 
5. Actually, V[A] is about 2276 times larger. Large variance values in a formulation using 
reliability bounds lead to the penalization of the compensator in spite of its actual improved 
performance. The reader should also notice that there is an offset of about 10% between 
E[A] and A(E[p]). Performance improvements in all the metrics can be seen in Figures 9 and 
10. 
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Figure 7. u(t) for the mean and variance base compensator 


B. Disturbance Rejection for a Flexible Beam 

The second example will focus on a reliability-based disturbance rejection solution for a 
flexible beam test article with both physical and modal parameter uncertainties. The system 
consists of a flexible thin aluminum blade, one-meter long, attached at its base to a hub 
motor. The hub motor is the control actuator for the system. At the tip of the beam, there 
is a reaction wheel that serves as a disturbance generator. The test article has nine sensors 
that may be used in any combination for either feedback or performance output monitoring. 
The finite element method is used to model this system by utilizing Euler-Bernoulli planar 
beam elements. A complete description of the flexible beam test article 20 is available. 

For this paper a SISO problem where the input u is the hub motor torque and the 
measured output y is the tip velocity, is studied. The tip reaction wheel disturbance is 
modeled by passing a Gaussian white noise process through a second-order linear low-pass 
filter, with parameters £ = 0.8 and a ; n = 2007T rad/s. The first five modes of the elastic 
structure are used to build a state space realization of the plant. This, in addition to 
the disturbance model leads to an open- loop system where x G M 12 , u G M. and y G M. 
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Table 2. Comparison of HSS and FSMSO for k^. 


V[p (i) ] for 
i = 1,...7 

Cost vector 

c 

Average 

error 

1 x 10“ 4 

C HSS = 10- 3 [0.000823, 0.0096, 0.5880, 0.0058, 0.0013] 
cfsmso = 10- 3 [0.000827, 0.0095, 0.5880, 0.0058, 0.0001] 

0.3% 

1 x 10~ 3 

c HSS = 10- 3 [0.00824, 0.0110, 0.6560, 0.0580, 0.0136] 
cfsmso = 10- 3 [0.00823, 0.0108, 0.7720, 0.0580, 0.0133] 

4% 

3 x 10~ 3 

c HSS = 10- 3 [5.1000, 0.0505, 0.8800, 0.1745, 0.0446] 
c fsmso = 10- 3 [0.0247, 0.0194, 1.9160, 0.1740, 0.0408] 

57% 

4 x 10“ 3 

c HSS = 10 _3 [32.400, 0.0013, 1.0880, 0.2330, 0.0626] 
cfsmso = 10~ 3 [0.0329, 0.0267, 2.8680, 0.2320, 0.0548] 

454% 


ME[p]) E[X(p)] 



Figure 8. CDF of A for a compensator with parameters kj. 


The uncertain parameters are the Young’s Modulus E (Pa), the density p (Kg/m 3 ) and the 
damping ratios of the retained vibration modes. This set leads to p = [. E , p, £ 2 , £ 3 , £ 4 , £s] T , 
whose components are assumed independent. The corresponding PDFs are given in Table 
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Figure 9. y(t) for a compensator with parameters kj. A zoom is shown below. 


3. The mean value of the parameters E[p] is set to coincide with the parameters in the 
finite element model. These mean values were chosen to match experimental data, while 
the supports of the distributions were set according to reasonable ranges of variation. The 
shapes of the PDFs were arbitrarily set. Performance requirements on stability and the 
output RMS are considered. Full-state feedback with a full-order observer determine the 
control structure. 

1. Nominal Compensator 

As before, a baseline compensator for the nominal plant is designed such that the RMS value 
is minimized. The resulting compensator, denoted as di, leads to y rm s = 0.011 m/s. The 
propagation of the uncertainty prescribed in Table 3 through the closed loop-system show 
that this compensator is robustly unstable with Fa(0) = 0.235. 
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Figure 10. u(t) for a compensator with parameters k 4 


2. Reliability -based Compensator 

For the sake of comparison, a reliability compensator will be synthesized using the hybrid 
method of reference . 13 In particular, a shapeable failure domain for the RMS requirement is 
assumed. This leads to the cost vector c = [r^(0), r yrms (e) + 7 yr ms) T , where e e [0, 0.05] and 
7 yrms = e - The selected control structure makes the feedback gain G, the observer gain L. 
and the RMS failure boundary e, the design variables. Recall that the separation principle 
does not hold. The resulting closed-loop dynamics is given by Equations (3) and (4). Notice 
that although the observer is deterministic, all the closed-loop poles are random. 

The synthesis approach with w = [20, 1] T leads to a compensator with parameters do, 
for which f A (0) = 0, e = 0.0139 m/s, f yrms (e) = 3.6 x 10 ~ 3 and c = [0,3.6 x 1CT 3 ] 2 ’. The 
probabilistic analysis of d 2 leads to Figures 11-12. Figure 11 shows that the random variable 
y rm s is moved toward zero, by virtue of the non- fixed failure boundary. In addition, Figure 
12 shows Bode magnitude plots of the disturbance to output transfer function, namely T zy . 
Notice that differences in the low-frequency portion of the diagram have a bigger impact 
on the RMS value. In addition, considerable variability in the closed-loop Bode magnitude 
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Table 3. Uncertainty model. 


E 

A e = 10 10 [5. 226, 7.839] 

ME) = B( 5,5) 

P 

A p = [2280, 3420] 

SM = B( 3,3) 

6 

A Sl = [0.08,0.12] 

f^) = B{2,2) 


A ?2 = [0.0252, 0.0378] 

£ 2 ( 6 ) = 5 ( 2 , 2 ) 

£3 

A ft = [0.02, 0.03] 

M&) = B( 2,2) 

£4 

A ?4 = [0.0304, 0.0456] 

MU) = B(2,2) 

£5 

A& = [0.02, 0.03] 

Mts) = B( 2,2) 


plot as well as a significant reduction in the damping of the first mode are attained near 20 
rad/s. It is interesting to notice that even though d 2 leads to a robustly stable closed-loop 
system in Equation (3), the full-state feedback subsystem A u and the full-order observer 
subsystem A 2j2 have a non-zero probability of instability. This indicates that the Separation 
Principle artificially reduces the design space once uncertainty is present. 



Figure 11. PDFs for the reliability based and the robustness based compensators. 
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Figure 12. Bode diagrams of T zy for d 2 - 


3. Mixed Compensator 

Lets take c = [f A (0), r yrms (0)] T , where the hybrid approach of reference 13 will be used for the 
stability metric and HSS for the RMS metric. The synthesis algorithm leads to d 3 , for which 
c = [0, 1.34 x 10~ 4 ] 7 . This compensator leads to the dashed line in Figure 11. Comparing 
both solutions, we see that while the PDF corresponding to the reliability-based compensator 
has less probability of exceeding the boundary value of e = 0.0139 m/s, the robustness-based 
compensator leads to a PDF which is much more concentrated toward the ideal value of zero. 
This clearly shows that the conceptual differences between the two formulations. Since there 
is no conservatism in the selection of the nominal plant, i.e., G(E[p]) is not the most difficult 
plant to control, y rms = 0.011 m/s does not necessarily bound the supports of the PDFs. This 
can be observed in Figure 11, where the support corresponding to d 3 contains y rms = 0.011 
m/s. 
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4- Mean and Variance based Compensator 

For this case we assume c = [6 A (0), T yrms (0)] T , where the FSMSO method will be used 
to estimate both indexes. In contrast to the satellite attitude example, the accuracy of 
the FSMSO method made it suitable for synthesis. Analytical sensitivities were used in 
the moments approximations. For instance, the sensitivities of the closed-loop poles that 
determine A are given by 

Avj = SjWj 

dsj = zJdAvj 

d 2 Sj = z 1 J d 2 A + (<9A — dsjl)(sjl — A) -1 <9A + dA(sjI — A) -1 (<9A — dsjl) Vj 

where z 3 is the jth eigenvector of A 1 , the right and left eigenvalues are normalized, i.e. 
zj 'Vj — 1, and and non- repeated poles are assumed. See the work of Burchett and Costello 21 
for a review on the subject. Derivatives of the output covariance are given by 

dy rms = jdiag G9QC T + 2<9CQC T } V2 

d 2 y rms = jdiag C<9 2 QC T + 4dCdQC T + 2<9 2 CQC T + 2dCdQdC T } V2 

where the derivatives of the state covariance are given by the solution to the set of Lyapunov 
equations 

A<9Q + dQA T + <9AQ + QdA T + dBSB T + BSdB T = 0 

A(9 2 Q + <9 2 QA T + 2dA8Q + 2dQdA T + c> 2 A9Q+ 
dQd 2 A T + 2dBSdB T + d 2 BSB T + BS<9 2 B T = 0 

The synthesis algorithm leads to d 4 , for which c = [5.33 x 10” 6 , 1.33 x 10~ 4 ] T . The resulting 
PDF for the RMS is indistinguishable from the one shown in Figure 11 even though the 
compensators are different. The CDFs of A for d 3 and d 4 are superimposed in Figure 
13. As before, the tendency of the formulation of making A as deterministic as possible is 
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apparent. Notice however, that the CPU time required for the analysis of d 3 is 421 s while 
the one for cU is 1.22 s. On the other hand, the analysis of d 4 via HSS takes 64 s. This 
exemplifies substantial savings in CPU time which result from using the FSMSO method. 
In this example, those savings justify the labor required to compute analytical derivatives. 
The reader must recall, however, that the same method led to inaccurate estimations for the 
satellite example. 
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Figure 13. CDFs of A for d 3 and d 4 . 
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VIII. Conclusions 

This paper proposes a control synthesis methodology for systems with probabilistic un- 
certainty. Synthesis is performed by solving a multi-objective optimization problem which 
combines requirements of stability and performance in time- and frequency-domains. In this 
study, reliability- and robustness- based formulations are proposed and several numerical 
methods for estimation are examined. In a reliability formulation, the probability of violat- 
ing design requirements is minimized while admissible domains are contracted toward regions 
with an improved performance. In a robustness-based formulation, a random metric that 
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measures the concentration of the random variable /pro cess about a target scalar/function is 
minimized. These two formulations lead to compensators with distinctive characteristics. In 
addition, metrics that bound the reliability metrics proposed, whose estimation only requires 
of means and variances, are also derived and used for control design. 

Some of the fundamental differences between the proposed strategy and conventional ro- 
bust control methods are: (i) unnecessary conservatism is eliminated since there is not need 
for convex supports/sets, (ii) the most likely plants are favored during synthesis allowing for 
probabilistic robust optimality, (iii) the trade off between robust stability and robust per- 
formance can be explored numerically, (iv) the uncertainty set, which could be unbounded, 
is closely related to parameters with clear physical meaning and (v) compensators with im- 
proved robust characteristics for a given control structure can be designed, e.g. one can 
search for a PID controller with best robust characteristics. Examples of attitude control of 
a simple satellite model and disturbance rejection of a flexible beam are used to elucidate 
the nature of the problem at hand and to demonstrate the methodology. 
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